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ABSTRACT 

We study the matter and velocity divergence power spectra in a f(R) gravity theory and their 
time evolution measured from several large-volume iV-body simulations with varying box 
sizes and resolution. We find that accurate prediction of the matter power spectrum in f(R) 
gravity places stronger requirements on the simulation than is the case with ACDM, because 
of the nonlinear nature of the fifth force. Linear perturbation theory is shown to be a poor ap- 
proximation for the f(R) models, except when the chameleon effect is very weak. We show 
that the relative differences from the fiducial ACDM model are much more pronounced in 
the nonlinear tail of the velocity divergence power spectrum than in the matter power spec- 
trum, which suggests that future surveys which target the collection of peculiar velocity data 
will open new opportunities to constrain modified gravity theories. A close investigation of 
the time evolution of the power spectra shows that there is a pattern in the evolution history, 
which can be explained by the properties of the chameleon-type fifth force in f(R) gravity. 
Varying the model parameter \fno\, which quantifies the strength of the departure from stan- 
dard gravity, mainly varies the epoch marking the onset of the fifth force, as a result of which 
the different f(R) models are in different stages of the same evolutionary path at any given 
time. 
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1 INTRODUCTION 

The origin of the observed accelerated e xpansion of the Universe 
jRiess etal.ll 19981 ; I Perlmut ter et al. I ll999h is one of the most chal- 
lenging questions in contemporary theoretical physics. Although 
the standard cold dark matter (CDM) model plus a cosmologi- 
cal constant can explain this observation very well, the so-called 
ACDM paradigm suffers from serious theoretical problems, as the 
vacuum energy density predicted by particle physics theory is many 
orders of magnitude larger than the cosmologically inferred value 
of the cosmological constant. This has motivated the proposal of 
alternative models to explain the accelerated expansion of the Uni- 
verse. 

So far, most of these models can b e divided into two classes: 
dark energy (see ICopelandetaLll2006l for a review), which in- 
volves one or more dynamical fields or new matter species that 
accelerate the expa nsion of the Universe, and modified gravity 
JCliftonetai1l2012h , which proposes that general relativity (GR) 
breaks down on cosmological scales and must be accompanied by 
certain modifications. Other models, such as the inhomogeneous 



universe model dBiswas & Notarill2008l) . have also been studied as 
alternatives to ACDM, but to a lesser extent. 

Unlike pressure-less matter, usually dark energy does not clus- 
ter _stronj*by_(which is the case for, e.g., the quintessence model 
of IWang et ai]|2000h and its effects are mainly to modify the cos- 
mic expansion history (there a re, however, exce ptions, such as the 
coupled quintessence model of lAmend ola 2000, in which the dark 
energy field does experience strong clustering). In these models, 
structure formation is different from that in ACDM only because 
the background expansion rate has been modified. In contrast, mod- 
ified gravity models often predict a different force law between 
matter particles, therefore changing structure formation directly. 
One can therefore in principle distinguis h between these po ssibili- 
ties using a combination of observables dJain & Zhan3l2008l) . 

Any modification to the force law in modified gravity theo- 
ries is highly constrained, becaus e GR has been confirmed to high 
accuracy by local tests (see, e.g.. iHovle et al. I l200ll ; Bertotti et al.1 



2003; Adelberger et al. 20031; iLvne et al.l 12004; Will 2006). If we 
consider the modification to the standard gravity as a new force, 
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the so-called fifth force, then the fifth force must either have a 
very weak strength or very short (sub-millimeter) range, in order 
to be consistent with local tests. Consequently, any viable modified 
gravity theory must have some mechanism to suppress (or screen) 
the fifth force at least in regions where local tests have been car- 
ried out. In the case that the new force is mediated by a scalar de- 
gree of freedom, there are several elegant exa mples of such screen- 
ing mechanisms, includin g the ch ameleon ( Khourv & WeltmarJ 
|2004 iMota & ShawH2007l). dil aton dBrax et al.ll201Cl>. svmrnetron 
fainterbichler & KhountolOT) and Ve inshtein dDvali et alj|200(j| : 
iNicolis et alj2009l ; lDeffavet et alj|2009h . 

In this work, we focus on one of the most well-st udied modi- 
fied gravity models, f{R) gravity dCarroll et al.ll2005l) . which em- 
ploys the chameleon mechanism to suppress the fifth force in high- 
density regions. In particular, to study the behaviour of the matter 
and velocity divergence power spectra in f(R) gravity, we perform 
a number of TV-b ody simulations for the f(R) model proposed 
in iHu & Sawic ki (2007) with various model parameters and sim- 
ulation box sizes (see, e .g. Jstarob inskv 2007c iLi & Barrowll2007l ; 
lApplebv & Battvdl200l for som e other viable f(R) c osmological 
models studied in the literature.). I Jennings et al.l (2012) used some 
of these simulations to study the form of redshift space distortions 
in f(R) models. Here we study the power spectra behind these dis- 
tortions in more detail. 

The nonlinear matter pow er spectrum in this p artic ular model 
has be en studied previously bv lOvaizu et al.l ( 120081) and lZhao et al.l 
d201ll) . Our study differs from these in several aspects. Firstly, the 
simulations used here have the largest volume up to date for models 
of this type and span a wide range of box sizes with good agreement 
between the results of each simulation. Secondly, we present the 
measurements of the velocity divergence power spectrum. Thirdly, 
we study the time evolution patterns for both the matter and veloc- 
ity divergence power spectra, and relate them to the property of the 
chameleon fifth force and the formation of structure in hierarchical 
cosmologies. 

The outline of the present work is as follows: in §|2]we briefly 
describe the general f(R) gravity models and explain how the 
chameleon mechanism works. In §[3] we give a short description 
of the simulation code and the technical specifications of the sim- 
ulations. §|4]contains the main results of this paper and we finally 
summarise and conclude in § [5] 

Throughout this paper we adopt the unit convention c = 1, 
where c is the speed of light. 



2 THE F(R) GRAVITY THEORY 

This section is devoted to a brief overview of the f(R) gravity the- 
ory and its properties. 



2.1 The f(R) gravity model 

The f(R) gravity model dCarroll et all 120051) is a straightfor- 
ward generalisation of GR: the Ricci scalar R in the Einstein- 
Hilbert action is replaced with an algebraic function f( R)(see e.g., 
ISotrriou & FaraonfeoiOt Ide Felice & Tsuiikawdlioia for recent 



reviews): 



S = 



M 2 

^l[R + f(R)]+C„ 



(1) 



in which Mpi is the Planck mass, M p 



8nG, G is Newton's 



grangian density for matter fields (photons, neutrinos, baryons and 
cold dark matter). By designing the functional form of f(R) one 
specifies the f(R) gravity model. 

Varying the action Eq. (0 with respect to the metric g^ v yields 
the modified Einstein equation 



gn V - V M V„/ji = 8ttGT™ , (2) 



in which G M „ = R^ — ^g^R is the Einstein tensor, fa = 
d//di?, V M the covariant derivative compatible to the metric g^, 
□ = V°V a and T™ is the energy momentum tensor for matter. 
One can consider Eq. (|2j as a fourth-order differential equation, or 
alternatively the standard second-order equation of GR with a new 
dynamical degree of freedom, fa, the equation of motion of which 
can be obtained by taking the trace of Eq. ([2} 

nf R = ±(R-f H R + 2f + 8-KGp m ), (3) 

where p m is the matter density. The new de gree of freedom fa is 
sometimes dubbed scalaron in the literature dZhao et aill201 ID . 

Assuming that the background Universe is described by the 
flat Friedmann-Robertson-Walker (FRW) metric, the line element 
in the perturbed Universe is written as 

ds 2 = a 2 ( V ) [(1 + 2$)d?i 2 - (1 - 2<t>)dx l dx^ , (4) 

in which r\ and x 1 are respectively the conformal time and comov- 
ing coordinates, $(77, x) and ^(77, x) are the Newtonian potential 
and perturbation to the spatial curvature, which are functions of 
both time 77 and space x; a denotes the scale factor of the Universe 
and a = 1 today. 

As we are mainly interested in the large-scale structures much 
smaller than the Hubble scale, and since the time variation of fa 
is very small in the models considered below, we shall work in the 
quasi-static limit by neglecting the time derivatives of fa. In this 
limit, the scalaron equation reduces to 



V 2 .fa = -\a [R(f R )- R + 8nG( P m- Pr, 



(5) 



in which V is the three dimensional gradient operator (to be distin- 
guished from the V introduced above), and the overbar means the 
background value of a quantity. Note that R can be expressed as a 
function of fa. 

Similarly, the Poisson equation which governs the Newtonian 
potential $ can be simplified to 



16ttG 



a 2 (p„ 



pm) 



1 2 

—a 

6 



R(f R )-R] 



(6) 



constant, g is the determinant of the metric g M „ and C m is the La- 



by neglecting terms involving time derivatives, and using Eq. (|5j to 
eliminate vfa. 

According to the above equations, there are two potential ef- 
fects of the scalaron on cosmology: (i) the background expansion 
of the Universe may be modified by the new terms in Eq. (|2j and (ii) 
the relationship between gravity and the matter density field is mod- 
ified, which can change the matter clustering and growth of density 
perturbations. Clearly, when \f R \ <C 1, we have R m —8nGp m 
from Eq. (|5} and so Eq. ((6) reduces to the normal Poisson equa- 
tion in GR; when \f R \ is large, we instead have \R — R\ <C 
8nG\pm — pm\ and so Eq. ® reduces to the normal Poisson equa- 
tion with G rescaled by 4/3. Note that this 4/3 is the maximum 
enhancement factor of gravity in f(R) models, independent of the 
specific functional form of f(R). The choice of f(R), however, 
is important because it governs when and on which scale the en- 
hancement factor changes from 1 to 4/3: scales much larger than 
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the range of the modification to Newtonian gravity mediated by the 
scalaron are unaffected and gravity is not enhanced there, while on 
much smaller scales the 4/3 enhancement is fully realised - this 
results in a scale-dependent modification of gravity and therefore a 
scale-dependent growth rate of structures. 



2.2 The chameleon mechanism 

The f(R) model would have been ruled out by local tests of 
gravity due to the factor-of-4/3 enhancement to the strength of 
Newtonian gravity. F ortunately, it is well known that, if f(R) is 
chosen appropriately (j Brookfield et al. 20od [Faulkner et alll2007l: 
Navarro & Van Acoleyenll2007t iLi & Barrow|[2007l: iHu & Sawickil 
2007; Brax et al. 2008) , the mode l can exploit the cham eleon mech- 
anism (Kh ourv & Weltmardl2004l ; lMota & Shawll2007r ) to suppress 
the enhancement and therefore pass the experimental constraints in 
high matter density regions such as our Solar system. 

The essence of the chameleon mechanism is as follows: the 
modifications to the Newtonian gravity can be considered as an ex- 
tra, or fifth force mediated by the scalaron. Because the scalaron it- 
self is massive, the force is of the Yukawa type and is suppressed by 
an exponential factor exp(— mr), in which m is the scalaron mass 
and r the distance between two test masses. In high matter density 
environments, m is very heavy and the suppression becomes very 
strong. In reality, this is equivalent to setting \/r\ <C 1 in high den- 
sity regions because of the exponential suppression, which leads to 
the GR limit as discussed above. 

As a result, the functional form of f(R) is crucial to determine 
whether the fifth force can be sufficiently suppressed in high den- 
sit y environments . In th is work we study the f(R) model proposed 
for which 



41M 3> M , and this simplifies the expression of the scalaron to 



sit y environments . In thi s 
bv lHu& Sawickil d2007l) . 



f(R) = -M 5 



Cl 



-R/M 2 



c 2 (--R/M 2 )" + l' 



(7) 



where M 2 = 87rGp m o/3 = Hgtlm, where H is the Hubble ex- 
pansion rate and Q m is the present-day fractional density of mat- 
ter. Hereafter a subscript o always mea ns the present day (a = 1) 
value of a quantity. It was shown by IHu & Sawickil ( 120071) that 
| /ho | < 0.1 is required to evade the Solar system constraints but 
the exact value depends on the behaviour of f R in galaxies as well. 

In the background cosmology, the scalaron fn always sits 
close to the minimum of the effective potential that governs its dy- 
namics, defined as 

V eB (f R )=~{R-f R R + 2f + 8wGp m ), (8) 

aroun d which it quickly oscillates with small amplitude dBrax et all 
|2012|) . Therefore we have 



R » 8nGp m - 2/ = 3M 2 ( a" 3 + — 

3C2 



(9) 



To match the ACDM model in background evolution, we need to 
set 



(-2 



(10) 



where fi m (f2A) is the present day fractional energy density of the 
dark matter (dark energy). 

By taking f2 A = 0.76 and Q, m = 0.243, we find that \R\ « 



ci ( M 



2 \ "+ 1 



-R 



(11) 



Therefore, two free parameters, n and ci/c|, completely specify 
the f(R) model. Indeed, the latter is related to the value of the 
scalaron today, f R o, as 



Cl 



3 1 + 4 



(12) 



In what follows we study three f(R) models with n = 1 and 
I/hoI = 10" 6 , 10" 5 , 10~ 4 , which will be referred to as F6, F5 
and F4 respectively. These choices of | f R o | are meant to cover the 
whole parameter space that is cosmological interesting: if \fno\ > 
10 -4 th en the f(R) model violates the cluster abundance con- 
straints dSchmidt et a"i]|2009h . and if \f R0 \ < 10~ 6 then the differ- 
ence from ACDM would be too small to be observable in practice 
(see our results presented below). 



3 A-BODY SIMULATIONS OF F(R) GRAVITY 

From Eqs. 10 [6) we have seen that, given the matter density field, 
we can solve the scalaron field f R from Eq. (|5} and plug it into the 
modified Poisson equation (fSJ to solve for $. Once $ is at hand, 
we can difference it to calculate the (modified) gravitational force 
which determines how the particles move subsequently. That is ex- 
actly what we need to do in A-body simulations to evolve the mat- 
ter distribution. 

The main challenge in A-body simulations of models such 
as f(R) gravity is to solve the scalaron equation (|5}, which is in 
general highly nonlinear. One way to achieve this is to use a mesh 
(or a set of meshes) on which f R could be solved. This implies 
that mesh-based Af-body codes are most convenient. On the other 
hand, tree-based codes are more difficult to apply here, as we do 
not have an analytical formula for the modified force law (such as 
r~ 2 in the Newtonian case) due to the complexities stemmed from 
the breakdown of the superposition principal, or the invalidity of 
Birkhoff theorem in modified gravity. 

Af-body simulations of f (R) gravity an d related theories hav e 



previously been performed by Ovaizull 20081): Ovaizu et al 



I Schmidt et al.l |2009h:lzhao et alj d201 ll):lLi & Zhad feOOS 



(2008) 



2010); 



Schmidtl d2009i):lLi & B arrow! d201 ll) : iBrax et al.l d201 ll) ; ILi & Hul 



d201 ll) ; fpavis et all d2012h . However, these simulations are mostly 
limited by either the box size or resolution, or both. For this work 
we ha ve run simulat ions using the recently developed ECOSMOG 
code dli et al.ll2012l) . ECOS MOG is a modi fication of the mesh- 
based A^-body code RAMSES dTevssieiteOOll) . which calculates the 
gravitational force by first solving the Poisson equation on meshes 
using a relaxation method to obtain the Newtonian potential and 
then differencing the potential. The code does not solve gravity 
by summing over the forces from nearby particles exp licitly, such 
as tree-based codes like GADGET dSpringel et al.l200ll) . Additional 
features of the ECOSMOG code include: 

(i) The adaptive mesh refinement (AMR), which refines a mesh 
cell (i.e., splits it into 8 children cells) if the number of particles 
in a cell exceeds a pre-defined number (the refinement criterion). 
This gives a higher force resolution in high matter density regions 



1 These values are used in the f(R) simulations extensively in the litera- 



ture, and are adopted in the simulations of this paper in order to compare 
with previous work. 
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where the chameleon effect is strong and the f(R) equation is more 
nonlinear. The refinement criterion is normally chosen as a number 
between 8 and 12, and in our simulations we use 9. We find that 
this refinement criterion works well in our case, namely, it gives the 
required force resolution without generating an overly large com- 
putational overhead. 

(ii) The multigrid relaxation algorithm that ensures quick con- 
vergence. The relaxation method finds the solution to an elliptical 
partial differential equation (PDE) on a mesh by iteratively updat- 
ing the initial guess until it converges, i.e., becomes close enough 
to the true solution. But the rate of converges slows down quickly 
after the first few iterations. To improve on this, one can coarsify 
the PDE, i.e., move it to a coarser mesh, solve it there and use the 
coarse solution to improve the solution on the original fine mesh. 
Unlike other codes, ECOSMOG does this on all the AMR meshes, 
greatly improving the convergence behaviour over the whole com- 
putational domain. 

(iii) The massive parallelisation which makes the computation 
very efficient. This is the key feature that enables us to run large 
simulations such as the ones employed in this study, which are 
beyond the reach of serial codes, like the ones developed by 
iLi & Zhaol j200l,l2010l) ; lLi & Barrowl d201 ll) . 



A convergence criterion is used to determine when the relaxation 
method has converged. In ECOSMOG, convergence is considered to 
be achieved when the residual of the PDE, i.e., the difference be- 
tween the two sides of the PDE, is smaller than a predefined param- 
eter e. We have checked that for e < 10 ~ 8 the solution to the PDE 
no longer changes significantly when e is further reduced, and our 
c hoices of e will be listed in Table Q] Further details can be found 
in lLi etaljj2012h . 

In this work we study the matter density and velocity diver- 
gence power spectra in the f(R) cosmology over a wide range of 
scales and redshifts using simulations. All our numerical experi- 
ments are described by the same set of cosmological parameters, 
i.e. the background cosmology for all models is the same. The 
values of cosmological parameters for our runs are the following: 
Q m = 0.24, Q A = 0.76, h = 0.73, n s = 0.958 and <t 8 = 0.77. 
The first two are the present day values of the dimensionless en- 
ergy density of the non-relativistic matter (including baryonic and 
dark) and dark energy, h is the dimensionless Hubble parameter to- 
day, n 3 is the scalar index of the primordial power spectrum and as 
is the linear rms density fluctuation measured in spheres of radius 
8/i _1 Mpc at z = 0. All models in each simulation share the same 
initial condition ccomputed at the initial time of Zj = 49 using the 
Zel'dovich approximation ( Zel' dovichll 1970h . Note that in general 
the modified gravity af fects the generation of the initial condition 
too jLi & Barrowll20T ]]). Here we use the same initial conditions 
for all models within a given simulation set because the differences 
in clustering between GR and all our f(R) models are negligible 
at the starting redshift. 

The fact that we use the same initial conditions for simulations 
in a given set is an advantage. Since the initial density fields for the 
GR and f(R) simulations have the same phases, any difference in 
the power spectra that we find at later times will be a direct conse- 
quence of the different dynamics between the two cosmologies. We 
give more details describing our numerical experiments in Table[T] 



4 THE POWER SPECTRA OF F(R) GRAVITY 

We start by introducing the dark matter density field, given by the 
expression 

p(x,t) = (p) (1 + 5), (13) 

where (/>(£)) is the ensemble average of the dark matter density 
at time t, and S(x, t) describes local deviations from homogene- 
ity. Structure formation is driven only by the spatially fluctuating 
part of the gravitational potential, cj>(x,t), induced by the density 
fluctuation field 8. In f(R) cosmologies, however, we expect that 
in regions where the fifth force is not screened by the chameleon 
mechanism we will have an additional boost to the standard grav- 
itational potential induced by the scalaron as described by Eq. $6^. 
Thus we expect that to some extent clustering will be enhanced in 
our f(R) models. A convenient measure of the strength of dark 
matter clustering is the power spectrum. For a Fourier repre senta- 
tion of a real space density field 

5 n = (2tt)- 3/2 J S(x ) e"* 2 d 3 x , (14) 

the power spectrum is defined as (assuming spatial isotropy) 

P SS (k) EE P(k) = (|%| 2 ) . (15) 

In addition to the measure of dark matter clustering we are also 
interested in the statistical measure of the cosmic peculiar velocity 
field. The irrotational velocity field, v(x), can be characterised, up 
to an additive bulk velocity, by a single scalar field such as the 
velocity divergence 

6(x) = iv ■ v(x) . (16) 

The 9 is called the expansion scalar (see, e.g., |Peeblesll980h . Divi- 
sion of the velocity divergence by the Hubble constant makes this 
quantity dimensionless. The Fourier transform of the real space ex- 
pansion scalar is 

% = (2tt)" 3/2 J 9{x)e-^d 3 x , (17) 

and similary we can define the power spectrum of the velocity di- 
vergence 

Pee(k) = (j%| 2 ) . (18) 

In linear perturbation theory the ratio of the density power spec- 
trum to the power spectrum of the velocity divergence scaled by 
the square of the growth rate should be unity. However on nonlin- 
ear and weakly nonlinear scales this ratio deviates from unity as in 
the nonlinear regime velocities grow more slowly than the linear 
pertu r bation theory predict ion (see, e.g. JCiecielag & Chodorowskil 
l2004l ; |jennings et alj|201lh . Thus we would expect that the poten- 
tial effects induced by the f(R) fifth force in density and diver- 
gence power spectra may differ from one another in the nonlinear 
and weakly nonlinear regimes. 

4.1 The structure of the density and velocity fields 

We start our analysis by comparing the visual impression of the 
density and velocity divergence fields obtained for our GR and 
f(R) simulations. To measure the fields sampled from the dis- 
tribution of dark matter ( DM) particle positions and vel ocities 
we use the DTFE code of ICautun & van de Weygaertl (t2oTTh (see 
also ISchaap & van de Weygaertl l2000l : Ivan de Weygaert & S chaap 
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Table 1. Some technical details of the simulations performed for this work. F6, F5 and F4 are respectively the labels of the f(R) models w ith |/roI = 
10 — 6 , 10 — 5 , 10~ 4 . Here fcjvyg denotes the Nyquist frequency, e is the residual for the Gauss-Seidel relaxation used in the code <Li et all2 012). and the two 
values of the convergence criterion are for the coarsest level and refinments respectively. We also list in the last column the number of realisations for each 
simulation. 



models 


^box 


no. of particles 


k Nyq [/l/Mpc] 


force resolution [h 1 kpc] 


convergence criterion 


realisations 


ACDM. F6, F5, F4 


1.5h -1 Gpc 


1024 3 


2.14 


22.9 


\e\ < io~ 12 /io- 8 


6 


ACDM, F6, F5, F4 


1.0fc -1 Gpc 


1024 3 


3.21 


15.26 


|e| < 10- 12 /10- 8 


1 


ACDM, F6, F5, F4 


500/i- 1 Mpc 


512 3 


3.21 


30.52 


|e| < 10- 12 /10- 8 


1 


ACDM, F6, F5, F4 


250h _1 Mpc 


512 3 


6.43 


7.63 


|e < 10- 12 /10- 8 


1 



GR 



F5 




F6 



F4 



5 + 1 



Figure 1. (Colour Online) All model comparison of z = density fields (p m /p m = 1 + 5) for the 250ft -1 Mpc box. Each panel show a very thin slice 
(~ 0.5/i -1 Mpc) through the DTFE density field. The top panels are results for GR (left) and F6 (right), and the bottom panels show the results for F5 (left) 
and F4 (right) respectively. 



2009). The Delaunay tessellation has the advantage that the veloc- 
ity divergence field computed in this way is volume averaged, as 
required by calculation, rather than mass averaged. The DTFE code 
also avoids the problem of empty cells, where no particles can be 
found, which can arise i n direct assignment methods to measure the 
velocity field (see, e.g.. |Pueblas & Scoccimarroll2009l) . Hence the 
velocity divergence power spectrum measured in this way is unbi- 



ased and has better noise properties than is the case for standard 
interpolation methods that deal with mass weighted velocities. 

In Fig. [T] we plot thin slices (0.5/i -1 Mpc) from the z — 
DM density field in the 250ft" 1 Mpc box. The panels show the GR 
model (top-left), the F6 (top-right) and the F5, F4 models (bottom 
row from left to right). We note the large-scale structures and the 
patterns of the cosmic web are the same in all panels. This is not 
unexpected as these simulations started from the same initial con- 
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Figure 2. (Colour Online) Comparison of the velocity divergence fields for the GR (left panels) and F4 (right panels) models. Each panel shows a thin slice 
from the 250/i — 1 Mpc box, and each row corresponds to a different cosmic time as labeled: a = 1 (top), a = 0.7 (middle) and a = 0.5 (bottom). 
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ditions, thus they share the same phases. However, the delicate ef- 
fects of the f(R) fifth force can be observed on small scales, where 
some of the density field features like filaments and clusters appear 
thicker in f(R) universes. This is accompanied by deeper density 
dips in voids. 

The picture becomes even more interesting when we look at 
the velocity divergence fields, which are plotted in Fig.|2] The slices 
shown correspond to the same regions plotted in Fig. [T] from the 
B = 250/i -1 Mpc box. Here we plot only the GR and F4 models 
(for which the effect of the fifth force is strongest) for three distinct 
epochs of cosmic evolution: a — 0.5 (the bottom row), a — 0.7 
(the middle row) and a — 1.0 (the top row). From the plot we 
observe that the differences between the GR and F4 models are 
very small at earlier stages of evolution (the a — 0.5 case), but 
they become very prominent as we move towards the present day 
(a = 1). We note in particular that the velocity fields around and in- 
side filaments and clusters are characterised by higher divergence. 
It is clear from this figure that, in the cosmic webs (clusters and 
filaments) of both the GR and and f(R) cosmology, the velocity 
divergence shows larger deviations between the two cosmologies 
compared to the density field. 

Note that the velocity divergence field is positive in voids be- 
cause matter flows from the void centres and the flow increases 
near the edge closer to higher density regions. Because the fifth 
force speeds up the matter flow, the velocity divergence in voids 
is larger in F4 than that in GR, and this is clearer at earlier times 
(note the difference in the colours of the two bottom panels). The 
trend is reversed around the clusters because matter flows inwards 
here, and again the magnitude of 9 is larger for F4 due to the fifth 
force, which explains why the filaments appear to be thicker for F4 
in Fig. [2] Inside clusters and filam ents, the velocity divergence be - 
comes positive again, as noted bv lPueblas & Scoccima rro (2009). 
This is because of the virialisation: after the halo is formed, the 
particles stop falling into the potential wells but circle around the 
halo center, which can in principle make the velocity divergence 
positive. 

In the following sections we will precisely quantify and anal- 
yse these effects by studying both the matter and the velocity diver- 
gence power spectra in all ACDM and f(R) cosmologies. 

4.2 Measurement of power spectra 

The matter power spectrum has been measured fro m all the simu- 
lation s listed in Table[Tjusing two codes: POWMES (Colo mbi et al.l 
2009) and our own code that uses fields obtained from the DTFE 
method, which rely on different algorithms. POWMES constructs the 
density field on a regular grid by direct particle assignment, while 
DTFE first samples the density and velocity divergence fields us- 
ing Delaunay tessellation and then interpolates onto a regular grid. 
POWMES attempts to correct for the impact of the scheme used to 
assign particles to the FFT grid, whilst we do not attempt any such 
correction with the code using the DTFE method. 

The grid that is used for the power spectra measurement is 
chosen to have the same resolution as the domain grid in the simu- 
lations. For example, for the Lbox = 1.5ft Gpc and 1.0ft Gpc 
simulations the FFT grid has 1024 3 cells. Only the results of the 
Lhox = 1.5ft Gpc simulations have error bars, which show the 
scatter amongst all six realisations. 

To show the accuracy of our simulations and the power spec- 
trum codes, we plot the measured matter powe r spectra Pss from 
the ACDM simulations against the HALOFIT jSmifh et ai]|2003h 
prediction in Fig. [3] The HALOF I T result is obtained using the pub- 



licly available CAMB code jLewis et alj|2000h . assuming the same 
cosmological parameters as used in the simulations, and is used 
here simply as a reference. 

The left panel of Fig. [3] plots the A-body results measured 
using POWMES, and the HALOFIT power spectrum is plotted as 
a black solid line. It can be seen that the two agree very well 
over a wide range of length scales. For example, the matter power 
spectrum from the Lbox = 1.5ft ~ Gpc simulation is accurate 
for 0.004 < fc/(ft/Mpc) < 1. This is of course as expected 
given that the RAMSES code has been tested in many ways, and 
it gives us some reassurance about the ECOSMOG gravity solver. 
In f(R) gravity simulations, however, this relationship should only 
be used as a rough guide with caution, as we shall explain below in 
Sect.l431 

The right panel of Fig. [3] shows the corresponding spectra 
measured from the DTFE-constructed density and velocity fields. 
The results also agree with HALOF I T very well, especially on large 
scales. Note that in this case Pss starts to deviate from HALOFIT 
at smaller k. This is due primarily to the lack of any correction 
in this estimate for the effects of the scheme used to interpo- 
late the smoothed d ensity field onto the FFT grid (see |jindl2005l : 
IColombi et alj2009h . 

4.3 Resolution issues in f(R) simulations 

The shape of the matter power spectrum is sensitive to changes in 
the cosmological model assumed and as such it may be a sensitive 
probe of the underlying theory of gravity or the properties of dark 
energy. Here, we are mostly interested in the shape of APss / Pss, 
where APss is the difference between the matter power spectra for 
f(R) gravity and ACDM, defined as 

apss_ p/rw 

Pss -P s L s CDM (k) L - Uyj 

To make the plots clearer, we have rebinned. The data points are at 
the centres of the bins, and the average value of P(k) and error bars 
are computed as follows: for the Lbox = 1.5/i _1 Gpc simulations, 
the P(k) value is the average over all points in a given bin in all six 
realisations and the error bar shows the scatter amongst all these 
points; for all other simulations, the same thing is done but only for 
points within a given bin in a single realisation. 

The results for models F4/F5/F6 at a = 1 are shown in Fig. [4] 
and different symbols are used to denote different simulation box 
sizes. We can see that the different symbols overlap with each other 
quite well, especially on large scales. On small scales, the large-box 
simulations predict slightly larger difference in Pss, but the differ- 
ence is, in general, quite small, at least at this particular cosmic 
time. In Fig.|4]we have also overplotted the results from linear per- 
turbation calculation and HALOFIT. The HALOFIT results are ob- 
tained from the linear power spectra in f (R) gravity in the fitting 
formulae obtained by dSmifh et alj|2003l) . We have checked that 
on large scales t he HALOFIT result agrees with third-order per- 
turbation theory dKovamaetai] r2009) quite well, and both show 
better agreement with simulations compared with the linear per- 
turbation theory. Indeed, linear theory breaks down on almost all 
scales where the f(R) model deviates from GR, especially in the 
F4/F5 cases where the nonlinearity is stronger. On the other hand, 
the HALOFIT gives a worse fit to simulation results in F6. This is 
because HALOFIT is calibrated using GR simulations and it does 
not capture the effect of the chameleon mechanism in F6, which 
suppresses the deviation from ACDM on small scales. 

The shape of APss /Pss for the three f(R) models look 
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Figure 3. (Colour Online) The power spectra of the ACDM model from different simulation boxes (symbols as explained in the legend; 'Bxxx' means that 
the box size is xxx/i _1 Mpc) compared to the HALOFIT (black solid curve). Left panel: the matter power spectra measured using POWMES. Right panel: the 
matter (upper symbols) and velocity divergence (lower symbols) power spectra measured from the DTFE-constracted density and velocity divergence fields. 



quite different from each other. In the F6 case, this increases 
all the way down to the smallest scales probed by the simula- 
tions; in F5, a small bump appears in between k — 1/i/Mpc and 
k = 2/Mpc, while after that APgg/Pgg increases on small scales; 
for F4, a single peak appears at k ~ 1/i/Mpc, and on smaller 
scales APgg/Pgg simply decreases. In addition, the amplitude of 
A Pss/Pss increases with | fno \ ■ These features were also seen 
by iLi & Zhacl 420091 l201Ch in their chameleon simulations, then 
IZhao et alj fconri " in their small f(R) simulations, and are con- 
firmed here by our larger simulations; Recently, similar features 
have also been fo und in simulatio ns of generalised dilaton and sym- 
metron models jBrax et al.ll2012lFI . This can lead to the conclusion 
that in general the dynamics with a fifth force employing a spe- 
cific spatial screening mechanism leaves a characteristic mark vis- 
ible in statistics of spatial clustering. A proper understanding of 
the features which appear in the power spectra of our models can 
only be achieved when we have a clear picture of the time evolu- 
tion of APgg I Pgg, and for this we plot the results of F4 and F5 at 
a = 0.3,0.5 and 0.7 in Fig. [5] 

Before going through the details of the time evolution of 
APss I Pgg, let us note that our previous observation, that the large- 
box simulations tend to overestimate APgg/Pss on small scales, 
becomes more prominent at earlier times. As an example, at a — 
0.5 we find that APgs/Pgg has a peak value of ~ 70% accord- 
ing to the Lbox = 1.5/i _1 Gpc simulation, while this value de- 
creases for high-resolution simulations and drops to ~ 40% for the 
250fo~ 1 Mpc boxes. 

One probable reason for this is the different force resolutions 
in these simulations. As the particles cluster to form structures, lo- 
cal over-densities grow and the chameleon effect starts to suppress 
the fifth force. If the force resolution is too low, the density field 



tends to be underestimatecjf] and the fifth force overestimated, re- 
sulting in more clustering of matter in these simulations compared 
to the ones which have higher resolution. 

Note that this resolution effect is a separate issue which hap- 
pens for the calculation of the chameleon fifth force: with the same 
background density, putting a particle at the centre of a sphere with 
radius 7? or spreading its mass uniformly in the sphere produce the 
same gravity at R, but the fifth forces at 7? would be quite different 
in these two configurations. This has nothing to do with the resolu- 
tion that is required by the usual gravity (Poisson equation) solver. 

The implications of this result are: 

(i) Comparing the ACDM power spectrum with HALOFIT as 
we did in Fig.[3]does not automatically provides a useful guide as 
to the scales down to which the f(R) simulations are reliable. The 
resolution effect on the fifth force solver is complicated; it depends 
on the model (e.g., compare F4 and F5 at a — 0.5) as well as on 
the redshift (e.g., compare F4 at a — 0.5 and F4 at a = 1.0). 

(ii) One has to be careful in choosing the right resolution to ob- 
tain accurate results in the f(R) simulations. The most straightfor- 
ward way would be to run ever higher-resolution simulations and 
see where the results start to disagree with each other. 

As an illustration, for F4 the Lbox = 1.5/i _1 Gpc simu- 
lations can be trusted down to k ~ 0.4ft/Mpc at a = 0.5, 
k ~ 1/i/Mpc at a = 1.0, while its qualitative predictions can 
be trusted down to even smaller scales. Evidently, such a box is 
not good enough to study small structure in the simulations, but 
is sufficient to study large-scale properties, such as redshift-space 
distortions (Jennings et al.ll2012h . 



2 Other authors in their studies of cosmologies employing simpler forms 
of fifth force have found s imilar features in the matter power spectra as 
we ha ve observed here (e.g. Hellwing & Juszkiewicz 2009; Keselm an et alj 
l2009h . 



3 The density in each cell is computed using a triangular-shaped cloud as- 
signment scheme, and using larger cells means that the mass of a particle 
will be more widely spread and thus the density peaks lower. 
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Figure 5. (Colour Online) The time evolution of APss /Pss for models F4 (upper panels) and F5 (lower panels). The left, middle and right panels are 
respectively the results at a = 0.3, 0.5 and 0.7 (the a = 1 results are shown in Fig. |4). 'Bxxxx' in the legend means that the simulation box size is 
xxxxh~ 1 Mpc, the horizontal dashed line is identically zero and the solid black curve is the linear perturbation prediction which agrees with the simulations 
better at earlier times. 



4.4 Time evolution of the spectra 

Let us now discuss on the time evolution of APss/Pss- From 
Figs. |4] and [3] we can see that as the Universe evolves, not only 
the magnitude but also the shape of APss/Pss changes. Taking 
the F4 model as an example: at a = 0.3 APss/Pss increases as 
one goes to smaller scales (at least until fc ~ KWMpc); at a = 0.5 
a peak develops at fc pca k ~ 2/i/Mpc while APss / Pss decreases 
for k > fcp Ca k; then the peak of APss/Pss shifts towards larger 
scales with fc pca k ~ 1/i/Mpc at a = 0.7 and fc pca k ~ 0.9/i/Mpc at 
a — 1. The F5 model behaves similarly but the peak only develops 
at a ~ 1. Likewise, the F6 model does not develop any peak in 
APss/Pss by a = 1. 

For the velocity divergence power spectrum Pee and the cross 
power spectrum Pse, we use those measured from the Lbox = 
250/i~ 1 Mpc simulations. The result is not sensitive to the mass res- 
olution, though this box size is a bit too sma ll, which means that the 
meas ured Pee can be ~ 5 — 10% higher (Pueb las"& Scoccimarrol 
2009). However, we are interested in the qualitative behaviour 
rather than accurate measurement of Pee, and this box enables us 
to go to smaller scales. 

Fig. [6] shows the behaviour of Pee on small scales (0.03 ^ 
fc/(/i/Mpc) ^ 3) in the different cosmologies. To make the curves 
clearer we have plotted fc 3 / (2iv 2 )Pee(k) instead of Pee- From this 
plot we can see that 

(i) there is a peak-dip-peak pattern on small scales, which agrees 
with what we have seen in the velocity divergence field in Fig. [2] 



(ii) not only the power spectrum is enhanced by the fifth force, 
but the peaks and dip also shift towards larger scales as the fifth 
force becomes stronger. This implies that the peak-dip-peak pattern 
observed in the plots develops as structures form, and could possi- 
bly be related to the characteristic scales of the structure formation 
at a given time, as we will show in the next subsection. 

Fig.|7]illustrates the time evolution of APee / Pee (upper pan- 
els) and APse/Pse (lower panels) for models F4 (left column) 
and F5 (right column). We can see some interesting features in 
these plots. Taking APee /Pee of the F4 model as an example, 
at early times (e.g., a = 0.2) this ratio increases with k un- 
til small scales; a dip then develops e.g., at fcdip ~ 4/i/Mpc at 
a = 0.3 and the dip shifts towards larger scales at late times. 
Meanwhile, a peak appears at k s < fcdip, while for k > fcdi p 
APee /Pee goes up again. Furthermore, a second and minor dip 
could develop immediately to the right of fcdip- The features and 
evolution pattern for APse/Pse are very similar to APee /Pee, 
both being significantly larger than APss /Pss djennings et al.l 
|2012|): this implies that local measurements of the velocity field 
dPatiri et al.l2012l;lHudson & Turnbulll2012l;|Pike & Hudson! 20051 ; 
lKosowskv & Bhattacharya 2009; lDavis et al 201 1 ) can indeed be 
a good probe of modified gravity. 

Of course, the complicated shape and evolution of APee / Pee 
must come from the chameleon fifth force, and we shall present an 
explanation of this below. 

The above evolution pattern suggests that the time evolution of 
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Figure 7. (Colour Online) The time evolution of APgg/ Pgg in F4 (upper left panel) and F5 (upper right model), along with the time evolution of APsg /Pgg 
for model F4 (lower left panel) and model F5 (lower right panel). In all panels at k = 0.1/i/Mpc the symbols are for times a = 0.2, 0.3, 0.5, 0.7 and 1.0 
from bottom to top. All results are measured from the DTFE-constructed density and velocity divergence fields for the L^ ox = 250h ~ 1 Mpc simulations. 



APss I Pss and APgg /Pgg for F5 (F6) is just a postponed version 
of that for F4. If this is the case, then there is a natural explanation 
for this, namely that the whole evolution pattern is an effect of the 
fifth force, which is suppressed until later times for smaller values 
of | /ho | • We will describe this in more detail in the next subsection. 
Note that this time-shifting effect can also be seen in the linear 
perturbation results for APss /Pss and APgg /Pgg shown in Fig. [8] 

4.5 A tale of two universes 

Suppose that there are two universes which are completely identical 
except for the underlying gravity. In universe I, standard general 
relativity applies and in universe II, the f(R) gravity applies. The 
evolution of structure is the same in the two universes up to quite 
late times, say z = 49 which is the starting time of our experiments, 
since the fifth force in f(R) gravity is vastly suppressed until then. 
The subsequent evolution can be divided into several stages: 
1) Stage (a), the fifth force in f(R) gravity begins to af- 
fect increasingly larger scales, starting from the smallest one. This 
speeds up the flow of matter, making the smallest structures form 
earlier through collapse in universe II than in universe I (such a 
boost in the rate of the formation of hiera rchical structures was re - 
ported for a general fifth-force models bv lLi & Zhao1(l2009ll2010h ; 



iHellwing et all 1201 (j) and lli & Barrowl 1201 lh ). As a result, one 
can see that both APss /Pss and APgg /Pgg increase towards small 
scales. 

2) Stage (b), the small structures have formed, and this corre- 
sponds to shell crossing in the spherical collapse model in universe 
II, while in universe I the same structure is still forming. The veloc- 
ity divergence inside the collapsed regions becomes less negative 
and then positive, during which process its magnitude is smaller 
than in the collapsing regions [cf. Fig.[2]- Although this is the effect 
for an individual structure, we would expect to see this statistically 
(and same for the discussions below), i.e., in the power spectrum, 
because structures form earlier in F4 in general, and as a result a dip 
starts to appear in APgg / Pgg on the scale roughly corresponding to 
the collapsed structure. During this stage, APss /Pss still follows 
the pattern of stage (a). 

3) Stage (c), the same small structure collapses and forms in 
universe I as well. Inside the structure, the deepening of the total 
gravitational potential (with the fifth force contributing) in universe 
II makes matter move faster than it does in universe I. This fact is 
reflected as a continued increase in APgg /Pgg on scales larger than 
where the dip first appears. Meanwhile, the dip develops into a val- 
ley and moves towards larger scales, because structures form hier- 
archically and larger ones form later. During this stage, APgg /Pgg 
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Figure 8. (Colour Online) Linear theory predictions for the time evolution of APgg / Pss (left panel) and APgg/Pgg (right panel) for F4 (solid) and F5 
(dashed) at different expansion factors, as indicated by the labels. 



also continues growing on scales much larger than the collapsing 
regions. A small peak could appear on scales immediately smaller 
than that corresponds to the valley, because the divergence field 
crosses zero, making Peg smaller for GR; this is probably why there 
is a second and smaller dip in APgg /Peg. 

The evolution of APss /Pss in stage (c) is more complicated. 
During the first part, substage (ci), the small scale clustering is con- 
tinuously boosted by the fifth force in universe II, which means 
that APgg /Pss keeps growing towards small scales. Meanwhile, a 
bump starts to appear on scales roughly corresponding to the col- 
lapsing regions, reflecting the enhanced and earlier formation of 
larger structures in universe II, probably as well as the fact that the 
growth of APss /Pss on scales smaller than that corresponding to 
the bump is slowed down by the increased vel ocity dispersion in 
the structures in universe II jLi & BarrowfeOl lb . 

This continues into substage (cii), during which the bump in 
APss / Pss develops into a peak and shifts towards larger scales (at 
the same pace as the valley in APgg/Pgg shifts leftwards). At the 
same time APss /Pss goes down on small scales, because of the 
higher velocity dispersion inside halos in universe II. 

The above evolution history of the shapes of APss / Pss and 
APgg/Pgg depends on the properties of the fifth force (e.g., it 
grows in time), and could in principle be a unique feature of the 
chameleon-type modified gravity theories. According to this pic- 
ture, the different f(R) models studied in this paper should fol- 
low the same evolutionary path, but as the fifth force becomes non- 
negligible in different eras depending on the value of |/ro|, at any 
given time the evolution is at different stages for the different mod- 
els. 

As an illustration, at a = 1.0 F4 has reached stage (cii), F5 
reached stage (ci) and F6 somewhere between stages (b) and (ci). 
At a = 0.7, F4 and F5 are in stages (cii) and (b) respectively. At 
a — 0.5, F4 has just left stage (ci) while F5 is still in stage (b). 

Note that the decrease of APss / Pss towards small scales in 
the F4 model does not reflect the fact that the fifth force is sup- 
pressed in small systems such as clusters and galaxies, although 
the latter is true. 

Based on these observations, we can give a rough esti- 
mate of the scales where the peak in APss /Pss and the dip in 
APgg / Pgg appear. We define the variance of linear density fluc- 



tuations smoothed by a Gaussian filter as 

a 2 (R,z) = f kApL{ ^ ex p(~k 2 R 2 )dlnk, (20) 
J 2% 2 

where Pl (k) is the linear power spectrum. The characteristic mass 
of halos M* is defined by matching the variance of the linear den- 
sity fluctuation to the threshold density for collapse, S c , 

a{R„z) = 5 c , (21) 

where M, = 47ri?Jp m /3. We expect that on small scales at k > 
ft* = R^ 1 , the power spectrum is significantly affected by col- 
lapsed objects. Using the critical density obtained by the spherical 
collapse model S c = 1.673 for the ACDM model and 8 C = 1.692 
dSchmidt et alj2009l) for F^ll the characteristic scales are obtained 
as ft, = 1.05/iMpc" 1 in ACDM and ft, = 0.72/iMpc _1 in F4 at 
a = 1. At a = 0.5, these scales are give by ft* = 3.39h,Mpc _1 
in ACDM and ft* = 2.2/iMpc -1 in F4. The characteristic scale 
ft, — RZ 1 is always smaller in f(R) gravity as the nonlinearity is 
stronger than ACDM. From the above arguments, we expect that 
the peak in APss /Pss and the dip in APgg /Pgg appear roughly at 
ft* (F4) < k < ft* (ACDM), as on these scales, collapsed objects 
are already formed in F4, but they are still collapsing in ACDM. 
From Figs. [4] [5] [7] we can see that these scales are roughly consis- 
tent with the scales where the peak and the dip appear. We should 
emphasize that these scales only give qualitative estimates of scales 
where collapsed objects are important and it is necessary to study 
the formation of halos in detail to make more precise predictions. 



5 SUMMARY AND CONCLUSIONS 

To summarise, in this paper we have studied the shape and evo- 
lution of the matter and velocity divergence power spectra in the 
f(R) gravity model, with the aid of a number of high-resolution 
W-body simulations. For this we have work ed with the f(R) La- 
grangian proposed bv lHu & Sawic kHl l2007l) , fixing one of the two 

4 Note that here we have assumed that S c for F4 is scale-independent as in 
GR: this is obtained by rescaling the Newton constant by 4/3 everywhere 
and the chameleon effect is neglected. 



12 Baojiu Li et al. 



7 
6 
0.5 
0.4 



Q_ 0.3 

< 

0.2 



6<35 

0.30 
0.25 
0.20 
Q. 0.15 
0.10 
0.05 




Q_ 

< 



-0.05 
0.12 



Q. 

< 0.04 



B1500 H 'h%4 

B1000 

B500 

B250 

haloffi 

linear 



_i i i 




B1500 

B1000 

B500 

B250 

halofit 

linear 




data for z=0 from 1000 Mpc/h box 



0.01 0.1 1 10 

k (h/Mpc) 

Figure 4. (Colour Online) The relative difference between the matter power 
spectra of the f(R) and ACDM simulations at z = 0. The results have 
been binned along the fc-axis as described in the text. 'Bxxxx' in the leg- 
end means that the simulation box size is xxxxh~ 1 Mpc, and the horizontal 
dashed line is identically zero. The top to bottom panels show respectively 
results for models F4, F5 and F6. The black dotted and solid curves are 
respectively the pr edictions using linear perturbation theory and HALOFIT 
temith et alj2003l) . 



free parameters (namely setting n = 1 in Eq. [7}. This leaves us 
with only one free parameter \fno\, which is the present-day value 
of fn in the cosmological background. The value of |/r| controls 
the strength of the chameleon mechanism: the smaller \fn\ is, the 
stronger the chameleon effect becomes and the weaker the devi- 
ations from general relativity. Because \fn\ increases with time 
overall, a larger value of |/ho| means that the fifth force becomes 
unscreened at an earlier time. 

We have run a series of iV-body simulations to study the for- 
mation of cosmic structures in selected f(R) models using the 
ECOSMOG code. To assess all possible resolution and finite box 



Figure 6. (Colour Online) The patterns of the small-scale tails of the veloc- 
ity divergence power spectrum per octave k [i /2ir 2 Pgg(k) for the different 
cosmologies measured at z = from the B = 1000?t — 1 Mpc simulation 
box. 



effects we make sure that our simulations cover a wide range of 
length and mass scales. On very large scales, the matter power spec- 
trum of f(R) gravity is found to be the same as that of the ACDM 
paradigm, since these scales are well beyond the range of the fifth 
force. On small scales, the matter power spectrum develops non- 
trivial shapes, depending on the value of |/ho| and time. We stress 
that linear perturbation theory is a bad approximation even on large 
scales, especially for the cases with |/ro| = 10~ 5 and 10 -6 , in 
which the chameleon effect is strong and the scalaron equation is 
highly nonlinear. This implies that one should be cautious about 
forecasts made for modified gravity theories based on linear per- 
turbation theory calculations. In general full nonlinear numerical 
simulations are needed. 

The most challenging part of the f(R) simulation (and mod- 
ified gravity simulation in general) is that the fifth force becomes 
weak in high density regions, where higher resolution is needed. 
We have seen in § I4.3l that if the mass and force resolution is not 
high enough, the amplitude of density peaks could be underesti- 
mated and the magnitude of the fifth force overestimated, causing 
significant errors in the simulations. 

The peculiar velocity field in the f(R) gravity is more affected 
by the presence of the fifth force than the density field. Indeed, 
the velocity divergence power spectrum of the f(R) gravity can 
differ from that of ACDM by twice as much as the difference in 
the matter power spectrum (~ 100% versus ~ 50% for F4 and 
~ 60% versus ~ 30% for F5). Furthermore, the shape and evo- 
lution pattern of the velocity divergence power spectrum, although 
also dependent on | /jjo | and time, can be very different from those 
of the matter power spectrum. The large effect of modified grav- 
ity on the velocity divergence power spectrum, especially on small 
scales, implies that the motion of particles and thus the dynamical 
state of the halos can be very different in modified gravity theories. 
Galaxy rotational curves, for example, can be modified and this ef- 
fect is important when interpreting observational data. The spin of 
dark matter halos, especially in l ow-density regions, can also be 
significantly faster dlxe et af]|2012h . 

The dependencies of the matter and velocity divergence power 
spectra on |/_ro | and time can be simplified if one understands them 
as the dependency on a single quantity - the fifth force. We have 
shown that the shapes of the power spectra for different | /ro | actu- 
ally evolve on the same path, but for models with smaller j fm \ the 
fifth force is suppressed until later times and the whole evolution is 
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delayed. For example, the F5 APee /Pee at a = 1 looks like the F4 
results at a = 0.7. 

We have presented an explanation of the shape and evolution 
of the power spectra based on this observation, according to which 
the valley and peaks in APee / Pee appear as a result of the fact that 
structures form earlier in f(R) gravity than they do in the A CDM 
mo del. This also explains th e observation in Zha o et al.ld201lh (and 
also lLi&Zhaol ho09. 2010) for other chameleon-type models) that 
APss I Pgs first increases as k increases and later develops a peak 
at the k corresponding to the size of dark matter halos. 

In this paper we have only focused on the theoretical aspects 
of the power spectra in f(R) gravity. The qualitative results here 
are expected to be quite general, and according to the theoretical 
picture similar things would be found in other modified gravity the- 
ories with screening mechanisms, such as the symme tron and dila- 
ton models or the ReBEL model (Nusser et al. 2005). Our analysis 
could be generalised to those models, and also connections could be 
made to observations by, for example, considering the weak lens- 
ing shear spectrum etc., to place constraints on the parameter | fno | . 
These issues will be left to future work. 
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